clear; clc;

k = 3;
c1 = 1e-9; c2 = 0; c3 = 0;

mu_list = linspace(-3*pi, 3*pi, 100);
eta_list = linspace(-3*pi, 3*pi, 100);

xi1_vec = zeros(length(eta_list)*length(mu_list), 1);
xi2_vec = zeros(length(eta_list)*length(mu_list), 1);

parfor idx = 1:length(mu_list)*length(eta_list)
    [i, j] = ind2sub([length(eta_list), length(mu_list)], idx);
    mu = mu_list(j);
    eta = eta_list(i);
    
    [xi1, xi2] = solve_xi_system(k, mu, eta, c1, c2, c3);
    xi1_vec(idx) = xi1;
    xi2_vec(idx) = xi2;
end

% 循环外reshape为二维矩阵
xi1_map = reshape(xi1_vec, length(eta_list), length(mu_list));
xi2_map = reshape(xi2_vec, length(eta_list), length(mu_list));


% 绘制xi1和xi2的二维热图作为参考
figure;
subplot(1,2,1);
imagesc(mu_list, eta_list, xi1_map);
set(gca,'YDir','normal');
colorbar;
title('\xi_1 distribution');
xlabel('\mu'); ylabel('\eta');

subplot(1,2,2);
imagesc(mu_list, eta_list, xi2_map);
set(gca,'YDir','normal');
colorbar;
title('\xi_2 distribution');
xlabel('\mu'); ylabel('\eta');

% --- 你之前定义的函数 ---
function [xi1, xi2] = solve_xi_system(k, mu, eta, c1, c2, c3)
    fun = @(x) [ ...
        2*x(1) - k*sin(x(1) + mu) + k*sin(mu) - c2 - c3;
        x(1) + x(2) - k*sin(x(2) + eta) + k*sin(eta) + c1 - c2
    ];
    
    x0 = [0; 0];
    options = optimoptions('fsolve', 'Display', 'off', 'TolFun', 1e-12);
    sol = fsolve(fun, x0, options);
    xi1 = sol(1);
    xi2 = sol(2);
end


